Import Data

library(ggplot2)
library(plyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(naniar)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(ggpubr)
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
## 
##     mutate
theme_set(theme_pubclean())


# reading in data
# set working directory 
full_data <- read.csv("~/Downloads/bimuno/full_data.csv")

# grouping data to average per cage

full_data_group <- group_by(full_data, cage.code) 

full_data_nested <-  summarise(full_data_group,
                               cagecode=first(cage.code),
                       model=first(model),
                       drug=first(drug),
                       Weight_W1=mean(Weight_W1),
                       Weight_W2=mean(Weight_W2),
                       Weight_W3=mean(Weight_W3),
                       Weight_W4=mean(Weight_W4),
                       weight_euth=mean(Euthanisation_Weight),
                       OF_Distance_totalcm=mean(OF_Distance_totalcm),
                       OF_Velocity_meancms=mean(OF_Velocity_meancms),
                       FST_preswim_strug=mean(FST_preswim_strug), 
                       FST_preswim_swim=mean(FST_preswim_swim),
                       FST_preswim_imm=mean(FST_preswim_imm),
                       FST_swim_strug=mean(FST_swim_strug),
                       FST_swim_swim=mean(FST_swim_swim),
                       FST_swim_imm=mean(FST_swim_imm),
                       EPM_close_freq=mean(EPM_close_freq),
                       EPM_close_secs=mean(EPM_close_secs),
                       EPM_open_freq=mean(EPM_open_freq),
                       EPM_open_secs=mean(EPM_open_secs),
                       EPM_open_latency=mean(EPM_open_latency),
                       EPM_fullopen_freq=mean(EPM_fullopen_freq),
                       EPM_fullopen_secs=mean(EPM_fullopen_secs),
                       EPM_fullopen_latency=mean(EPM_fullopen_latency)
                       )

# remove NA columns
full_data_nested <- full_data_nested[rowSums(is.na(full_data_nested)) != ncol(full_data_nested), ]



#  Variable names: 
#  [1] "Cage.Code"            "animal.code"          "model"                "drug"                 
# "Weight_W1"            "Weight_W2"           "Weight_W3"            "Weight_W4"    "Euthanisation_Weight"   "OF_Distance_totalcm"  "OF_Velocity_meancms" 
#  [7] "FST_preswim_strug"    "FST_preswim_swim"     "FST_preswim_imm"      "FST_swim_strug"       "FST_swim_swim"        "FST_swim_imm"        
# [13] "EPM_close_freq"       "EPM_close_secs"       "EPM_open_freq"        "EPM_open_secs"        "EPM_open_latency"     "EPM_fullopen_freq"   
# [19] "EPM_fullopen_secs"    "EPM_fullopen_latency"

Testing Assumptions

###This experiment has two Factors (i.e., types of manipulations:  Phenotype and Invasiveness), and the experiment has data for all 4 possible combinations of these two Factors.  We should try to analyze it as a 2-way ANOVA.

###Plot the data:
full_data_nested$group <- factor(paste(full_data_nested$model,full_data_nested$drug))

boxplot(FST_swim_imm ~ group, data=full_data_nested,cex.axis = 1.2)
stripchart(FST_swim_imm ~ group, data=full_data_nested,
vertical = TRUE, method = "jitter",
pch = 21, col = "maroon", bg = "bisque",
add = TRUE)
mtext("immobility",2,line=2.5,cex=1.5)

###Model the data
p <- aov(FST_swim_imm ~ group, data=full_data_nested)

###Check the assumptions:
plot(p)

#### R identified case 13, 20 and 20 as outliers and have heterogeneity of residuals.. - all from Gin - Red group. 

leveneTest(full_data_nested$FST_swim_imm ~ full_data_nested$model * full_data_nested$drug)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.1862 0.3328
##       28
# Levene's Test for Homogeneity of Variance (center = median)
#       Df F value Pr(>F)
# group  3  1.1862 0.3328
#       28 

# so no significant differences in equal variances - therefore heterogeneity

shapiro.test(p$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  p$residuals
## W = 0.94601, p-value = 0.1109
#   Shapiro-Wilk normality test
# 
# data:  p$residuals
# W = 0.94601, p-value = 0.1109

# both not significant so no vioaltion of assumptions

hist(resid(p))

# redisuals are somewhat normally distributed

Results for FST Immobility:

All error bars throughout document are SD

## fst immobility
anova(lm(FST_swim_imm ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: FST_swim_imm
##            Df Sum Sq Mean Sq F value   Pr(>F)   
## model       1  10786 10786.1  8.0936 0.008213 **
## drug        1   4453  4453.3  3.3416 0.078221 . 
## model:drug  1    361   361.1  0.2710 0.606765   
## Residuals  28  37315  1332.7                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_fst_imm <- ddply(full_data_nested, c("model", "drug"), summarise,
               N    = length(FST_swim_imm),
               mean = mean(FST_swim_imm),
               sd   = sd(FST_swim_imm),
               se   = sd / sqrt(N)
)

# p1 <- ggplot(sum_data_fst_imm, aes(x=model, y=mean, fill=drug)) + 
#   geom_bar(position=position_dodge(), stat="identity") +
#   geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
#                 width=.2,                    # Width of the error bars
#                 position=position_dodge(.9))+
#                   labs(x = "", fill= "Drug")

# full_data_nested$dist_cat_n[full_data_nested$model == "Gin"] <- 1
# full_data_nested$dist_cat_n[full_data_nested$model == "Whisky"] <- 2
# 
# full_data_nested$scat_adj[full_data_nested$drug == "Blue"] <- -0.2
# full_data_nested$scat_adj[full_data_nested$drug == "Red"] <- 0.2
# 
# min.mean.sd.max <- function(x) {
#   r <- c(min(x), mean(x) - sd(x), mean(x), mean(x) + sd(x), max(x))
#   names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
#   r
# }


# p1 <- ggplot(full_data_nested, aes(x=model, y=FST_swim_imm, fill=drug))
# p1 +  geom_boxplot(
#   aes(color = drug), width = 0.5, size = 0.4,
#   position = position_dodge(0.8)
#   ) +
#   geom_dotplot(
#     aes(fill = drug, color = drug),
#     #trim = FALSE,
#     binaxis='y', stackdir='center', dotsize = 0.8,
#     position = position_dodge(0.8)
#   )+
#   scale_fill_manual(values = c("#0000FF", "#FF0000"))+
#   scale_color_manual(values = c("#0000FF", "#FF0000"))

#install.packages("Hmisc")
#library(Hmisc)
p_test <- ggplot(full_data_nested, aes(x=model, y=FST_swim_imm, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_test + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("FST Immobility (secs)")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

# 
# p1 <- ggplot(full_data_nested, aes(x=model, y=FST_swim_imm, fill=drug)) + 
#   geom_boxplot(outlier.size=0)+
#   labs(x = "", fill= "Drug")+ 
#   geom_jitter(aes(dist_cat_n + scat_adj,FST_swim_imm),
#         # position=position_jitter(width=0.1,height=0),
#         # alpha=0.6,
#         # size=3,
#         show.legend = FALSE)+
#     ylab("FST Immobility (secs)")+    stat_summary(fun.y=mean, geom="point", size=2, position=0.2) +
#   stat_summary(fun.data = mean_se, geom = "errorbar", position=-0.2)
# 
# p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))


# effect of model - expect to see difference between FSL and FRL in immobility time 
# no effect of drug - which is different to the non-averaged data

Results for FST Swimming:

## fst swimming
anova(lm(FST_swim_swim ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: FST_swim_swim
##            Df  Sum Sq Mean Sq F value  Pr(>F)  
## model       1  2236.1 2236.13  3.9288 0.05736 .
## drug        1  2153.3 2153.32  3.7833 0.06187 .
## model:drug  1   164.3  164.26  0.2886 0.59537  
## Residuals  28 15936.7  569.17                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_fst_swim <- ddply(full_data_nested, c("model", "drug"), summarise,
                          N    = length(FST_swim_swim),
                          mean = mean(FST_swim_swim),
                          sd   = sd(FST_swim_swim),
                          se   = sd / sqrt(N)
)

p1<- ggplot(sum_data_fst_swim, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")



p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:plyr':
## 
##     is.discrete, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
p_swim <- ggplot(full_data_nested, aes(x=model, y=FST_swim_swim, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_swim + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("FST Swimming (secs)")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

# no effect

Results for FST Struggling

## fst struggling

anova(lm(FST_swim_strug ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: FST_swim_strug
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## model       1  22845 22844.5 18.0540 0.0002148 ***
## drug        1    413   413.3  0.3266 0.5722185    
## model:drug  1   1012  1012.5  0.8002 0.3786659    
## Residuals  28  35430  1265.3                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_fst_strug <- ddply(full_data_nested, c("model", "drug"), summarise,
                           N    = length(FST_swim_strug),
                           mean = mean(FST_swim_strug),
                           sd   = sd(FST_swim_strug),
                           se   = sd / sqrt(N)
)

p1<-ggplot(sum_data_fst_strug, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

p_strug <- ggplot(full_data_nested, aes(x=model, y=FST_swim_strug, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_strug + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("FST Struggling (secs)")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

#model differences - no effect of drug
# attempt to create stacked bar chart of time spent in FST

# Activity <- c("Immobility", "Swimming", "Struggling")
# Model <- c("Gin Blue", "Gin Red", "Whisky Blue", "Whisky Red")
# 
# fst_stacked <- data.frame(full_data_nested$model, full_data_nested$drug, full_data_nested$FST_swim_imm, full_data_nested$FST_swim_swim, full_data_nested$FST_swim_strug)
# 
# fst_stacked$full_data_nested.model<- as.factor(fst_stacked$full_data_nested.model)
# fst_stacked$full_data_nested.drug<- as.factor(fst_stacked$full_data_nested.drug)
# fst_stacked$full_data_nested.FST_swim_imm<- as.numeric(fst_stacked$full_data_nested.FST_swim_imm)
# fst_stacked$full_data_nested.FST_swim_swim <- as.numeric(fst_stacked$full_data_nested.FST_swim_swim)
# fst_stacked$full_data_nested.FST_swim_strug <- as.numeric(fst_stacked$full_data_nested.FST_swim_strug)

# barplot(fst_stacked, main="Distribution of Amount of Time Spent in FST",
#   xlab="Groups", col=c("darkblue","red"),
#   legend = rownames(counts))

# library(lattice) 
# barchart(fst_stacked, scales = list(x = "free"), 
#           auto.key = list(title = "Time Spent in FST"), horizontal=FALSE)

Results for Open Field Total Distance:

#open field total distance

anova(lm(OF_Distance_totalcm ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: OF_Distance_totalcm
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## model       1 1786232 1786232 29.8951 7.736e-06 ***
## drug        1     202     202  0.0034   0.95408    
## model:drug  1  210284  210284  3.5194   0.07112 .  
## Residuals  28 1673000   59750                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_OF_Distance_totalcm <- ddply(full_data_nested, c("model", "drug"), summarise,
                            N    = length(OF_Distance_totalcm),
                            mean = mean(OF_Distance_totalcm),
                            sd   = sd(OF_Distance_totalcm),
                            se   = sd / sqrt(N)
)

p1 <- ggplot(sum_OF_Distance_totalcm, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

p_OF <- ggplot(full_data_nested, aes(x=model, y=OF_Distance_totalcm, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_OF + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("OF Total Distance Travelled (cm)")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

# model differences

Results for Open Field Speed/Velocity (m/s):

#open field speed velocity m/s

anova(lm(OF_Velocity_meancms ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: OF_Velocity_meancms
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## model       1 22.0110 22.0110 31.0170 5.867e-06 ***
## drug        1  0.0029  0.0029  0.0041   0.94919    
## model:drug  1  3.3168  3.3168  4.6739   0.03932 *  
## Residuals  28 19.8700  0.7096                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_of_velo <- ddply(full_data_nested, c("model", "drug"), summarise,
                            N    = length(OF_Velocity_meancms),
                            mean = mean(OF_Velocity_meancms),
                            sd   = sd(OF_Velocity_meancms),
                            se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_of_velo, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

# model differences - and model x drug interaction

Results for EPM Frequency of Full Body Entries to Open Arms:

#EPM full body to enter open arms frequency


anova(lm(EPM_fullopen_freq ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_fullopen_freq
##            Df  Sum Sq Mean Sq F value   Pr(>F)   
## model       1  34.031  34.031  8.5221 0.006852 **
## drug        1   0.125   0.125  0.0313 0.860841   
## model:drug  1   2.000   2.000  0.5008 0.484981   
## Residuals  28 111.812   3.993                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_epm_fullopen_freq <- ddply(full_data_nested, c("model", "drug"), summarise,
                          N    = length(EPM_fullopen_freq),
                          mean = mean(EPM_fullopen_freq),
                          sd   = sd(EPM_fullopen_freq),
                          se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_epm_fullopen_freq, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

p_EPM_full <- ggplot(full_data_nested, aes(x=model, y=EPM_fullopen_freq, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_EPM_full + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("EPM Frequency of Full Entries to Open Arms")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

# only model differences

Results for EPM Frequency of Entries to Closed Arms:

# EPM closed arms frequency to enter

anova(lm(EPM_close_freq ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_close_freq
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## model       1 118.195 118.195 34.1733 2.773e-06 ***
## drug        1  18.758  18.758  5.4234   0.02731 *  
## model:drug  1   2.820   2.820  0.8154   0.37422    
## Residuals  28  96.844   3.459                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_epm_closed_freq <- ddply(full_data_nested, c("model", "drug"), summarise,
                                    N    = length(EPM_close_freq),
                                    mean = mean(EPM_close_freq),
                                    sd   = sd(EPM_close_freq),
                                    se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_epm_closed_freq, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

#model and drug differences

Proportion time spent in open arms:

## make new variable time spent open / time spent open + closed arms
full_data_nested <- mutate(full_data_nested, EPM_proportion = (EPM_open_secs/(EPM_open_secs + EPM_close_secs)))


full_data_nested <- mutate(full_data_nested, EPM_proportion_total = (EPM_open_secs/300))

anova(lm(EPM_proportion ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_proportion
##            Df  Sum Sq   Mean Sq F value Pr(>F)
## model       1 0.01498 0.0149846  0.7037 0.4087
## drug        1 0.01367 0.0136674  0.6418 0.4298
## model:drug  1 0.00773 0.0077304  0.3630 0.5517
## Residuals  28 0.59626 0.0212950
sum_EPM_proportion <- ddply(full_data_nested, c("model", "drug"), summarise,
                                    N    = length(EPM_proportion),
                                    mean = mean(EPM_proportion),
                                    sd   = sd(EPM_proportion),
                                    se   = sd / sqrt(N)
)

p1 <- ggplot(sum_EPM_proportion, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

p_EPM_prop <- ggplot(full_data_nested, aes(x=model, y=EPM_proportion_total, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_EPM_prop + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("EPM Proportion of Time Spent in Open Arm (%)")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

## no sig differences


anova(lm(EPM_proportion_total ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_proportion_total
##            Df  Sum Sq   Mean Sq F value Pr(>F)
## model       1 0.01598 0.0159758  1.1460 0.2935
## drug        1 0.00850 0.0085043  0.6100 0.4413
## model:drug  1 0.00278 0.0027813  0.1995 0.6585
## Residuals  28 0.39033 0.0139405
sum_EPM_proportion_total <- ddply(full_data_nested, c("model", "drug"), summarise,
                                    N    = length(EPM_proportion_total),
                                    mean = mean(EPM_proportion_total),
                                    sd   = sd(EPM_proportion_total),
                                    se   = sd / sqrt(N)
)

p1 <- ggplot(sum_EPM_proportion_total, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

Results for EPM Frequency of Entries to Open Arms:

# EPM open arms frequency to enter

anova(lm(EPM_open_freq ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_open_freq
##            Df  Sum Sq Mean Sq F value  Pr(>F)  
## model       1  13.133 13.1328  3.4477 0.07389 .
## drug        1   2.258  2.2578  0.5927 0.44781  
## model:drug  1   6.570  6.5703  1.7249 0.19973  
## Residuals  28 106.656  3.8092                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_epm_open_freq <- ddply(full_data_nested, c("model", "drug"), summarise,
                                  N    = length(EPM_open_freq),
                                  mean = mean(EPM_open_freq),
                                  sd   = sd(EPM_open_freq),
                                  se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_epm_open_freq, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

# no significant differences

Results for EPM Time Spent with Full Body in Open Arms:

# EPM time spent full body in open arms in secs

anova(lm(EPM_fullopen_secs ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_fullopen_secs
##            Df Sum Sq Mean Sq F value Pr(>F)
## model       1   1770  1770.1  1.3293 0.2587
## drug        1   1540  1540.1  1.1566 0.2914
## model:drug  1    128   128.0  0.0961 0.7588
## Residuals  28  37286  1331.6
sum_data_epm_fullopen_secs <- ddply(full_data_nested, c("model", "drug"), summarise,
                                N    = length(EPM_fullopen_secs),
                                mean = mean(EPM_fullopen_secs),
                                sd   = sd(EPM_fullopen_secs),
                                se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_epm_fullopen_secs, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

p_EPM_total_open <- ggplot(full_data_nested, aes(x=model, y=EPM_fullopen_secs, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_EPM_total_open + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("EPM Time spent in Open Arms")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

# no significant differences

Results for EPM Time Spent with in Open Arms:

# EPM time spent in open arms in secs

anova(lm(EPM_open_secs ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_open_secs
##            Df Sum Sq Mean Sq F value Pr(>F)
## model       1   1438 1437.82  1.1460 0.2935
## drug        1    765  765.38  0.6100 0.4413
## model:drug  1    250  250.32  0.1995 0.6585
## Residuals  28  35130 1254.64
sum_data_epm_open_secs <- ddply(full_data_nested, c("model", "drug"), summarise,
                                    N    = length(EPM_open_secs),
                                    mean = mean(EPM_open_secs),
                                    sd   = sd(EPM_open_secs),
                                    se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_epm_open_secs, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

# no significant differences

Results for EPM Time Spent in Closed Arms:

# EPM closed arms time spent in seconds

anova(lm(EPM_close_secs ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_close_secs
##            Df Sum Sq Mean Sq F value Pr(>F)
## model       1     75   75.03  0.0551 0.8162
## drug        1    800  800.00  0.5871 0.4499
## model:drug  1   1188 1188.28  0.8721 0.3584
## Residuals  28  38152 1362.56
sum_data_epm_closed_secs <- ddply(full_data_nested, c("model", "drug"), summarise,
                                N    = length(EPM_close_secs),
                                mean = mean(EPM_close_secs),
                                sd   = sd(EPM_close_secs),
                                se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_epm_closed_secs, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

# no significant group differences

Results for EPM Latency to Enter Open Arms:

# EPM latency to enter open arms 

full_data_nested$EPM_open_latency <- as.numeric(full_data_nested$EPM_open_latency)

anova(lm(EPM_open_latency ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_open_latency
##            Df  Sum Sq Mean Sq F value Pr(>F)
## model       1   300.1  300.12  0.3157 0.5786
## drug        1   225.8  225.78  0.2375 0.6298
## model:drug  1  1300.5 1300.50  1.3682 0.2520
## Residuals  28 26614.8  950.53
sum_data_epm_open_latency <- ddply(full_data_nested, c("model", "drug"), summarise,
                                  N    = length(EPM_open_latency),
                                  mean = mean(EPM_open_latency),
                                  sd   = sd(EPM_open_latency),
                                  se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_epm_open_latency, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

# no significant differences

Results for EPM Latency to enter Open Arms with Full Body

# EPM latency to enter open arms with full body


anova(lm(EPM_fullopen_latency ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: EPM_fullopen_latency
##            Df  Sum Sq Mean Sq F value Pr(>F)
## model       1   182.9  182.88  0.2683 0.6085
## drug        1   130.0  130.01  0.1907 0.6656
## model:drug  1  1098.6 1098.63  1.6119 0.2147
## Residuals  28 19084.1  681.57
sum_data_epm_fullopen_latency <- ddply(full_data_nested, c("model", "drug"), summarise,
                                   N    = length(EPM_fullopen_latency),
                                   mean = mean(EPM_fullopen_latency),
                                   sd   = sd(EPM_fullopen_latency),
                                   se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_epm_fullopen_latency, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

## no significant differences

Weight:

anova(lm(Weight_W1 ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: Weight_W1
##            Df  Sum Sq Mean Sq F value   Pr(>F)   
## model       1  7595.3  7595.3 11.0173 0.002514 **
## drug        1  5671.1  5671.1  8.2262 0.007763 **
## model:drug  1  1785.0  1785.0  2.5893 0.118808   
## Residuals  28 19303.1   689.4                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_Weight_W1 <- ddply(full_data_nested, c("model", "drug"), summarise,
                                   N    = length(Weight_W1),
                                   mean = mean(Weight_W1),
                                   sd   = sd(Weight_W1),
                                   se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_Weight_W1, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

p_weight1 <- ggplot(full_data_nested, aes(x=model, y=Weight_W1, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight1 + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("Weight Week 1 (g)")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

Week 2

anova(lm(Weight_W2 ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: Weight_W2
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## model       1 18408.0 18408.0 30.8383 6.129e-06 ***
## drug        1  5791.6  5791.6  9.7024   0.00422 ** 
## model:drug  1  2476.3  2476.3  4.1485   0.05123 .  
## Residuals  28 16713.8   596.9                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_Weight_W2 <- ddply(full_data_nested, c("model", "drug"), summarise,
                                   N    = length(Weight_W2),
                                   mean = mean(Weight_W2),
                                   sd   = sd(Weight_W2),
                                   se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_Weight_W2, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

p_weight2 <- ggplot(full_data_nested, aes(x=model, y=Weight_W2, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight2 + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("Weight Week 2 (g)")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

Week 3

anova(lm(Weight_W3 ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: Weight_W3
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## model       1 23247.1 23247.1 33.4939 3.248e-06 ***
## drug        1  6146.6  6146.6  8.8560  0.005962 ** 
## model:drug  1  2269.7  2269.7  3.2701  0.081304 .  
## Residuals  28 19433.9   694.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_Weight_W3 <- ddply(full_data_nested, c("model", "drug"), summarise,
                                   N    = length(Weight_W3),
                                   mean = mean(Weight_W3),
                                   sd   = sd(Weight_W3),
                                   se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_Weight_W3, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

p_weight3 <- ggplot(full_data_nested, aes(x=model, y=Weight_W3, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight3 + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("Weight Week 3 (g)")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

Week 4

anova(lm(Weight_W4 ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: Weight_W4
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## model       1 26363.8 26363.8 41.0895 6.115e-07 ***
## drug        1  6626.9  6626.9 10.3284  0.003288 ** 
## model:drug  1  2354.7  2354.7  3.6699  0.065666 .  
## Residuals  28 17965.3   641.6                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sum_data_Weight_W4 <- ddply(full_data_nested, c("model", "drug"), summarise,
                                   N    = length(Weight_W4),
                                   mean = mean(Weight_W4),
                                   sd   = sd(Weight_W4),
                                   se   = sd / sqrt(N)
)

p1 <- ggplot(sum_data_Weight_W4, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

p_weight4 <- ggplot(full_data_nested, aes(x=model, y=Weight_W4, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight4 + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("Weight Week 4 (g)")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

Weight at Euthanisation

 anova(lm(weight_euth ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: weight_euth
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## model       1 25679.4 25679.4 40.0749 7.557e-07 ***
## drug        1  4548.2  4548.2  7.0978   0.01266 *  
## model:drug  1  3270.4  3270.4  5.1037   0.03185 *  
## Residuals  28 17942.0   640.8                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(full_data_nested$model, full_data_nested$drug, full_data_nested$weight_euth)

interaction.plot(full_data_nested$drug, full_data_nested$model, full_data_nested$weight_euth)

sum_data_weight_euth <- ddply(full_data_nested, c("model", "drug"), summarise,
                                   N    = length(weight_euth),
                                   mean = mean(weight_euth),
                                   sd   = sd(weight_euth),
                                   se   = sd / sqrt(N)
)

p1<- ggplot(sum_data_weight_euth, aes(x=model, y=mean, fill=drug)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9))+
                  labs(x = "", fill= "Intervention")

p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))

# p1 <- ggplot(full_data_nested, aes(x=model, y=weight_euth, fill=drug)) + 
#   geom_boxplot(outlier.size=0)+
#   labs(x = "", fill= "Drug")+ 
#   geom_jitter(aes(dist_cat_n + scat_adj,weight_euth),
#         # position=position_jitter(width=0.1,height=0),
#         # alpha=0.6,
#         # size=3,
#         show.legend = FALSE)+
#     ylab("Weight (g) at Euthanisation")
# 
# 
# p1+scale_fill_manual(values=c("#0000FF", "#FF0000"))



p_weight_euth <- ggplot(full_data_nested, aes(x=model, y=weight_euth, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight_euth + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("Weight (g) at Euthanisation")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

Weight Increase over 4 Weeks

# graph for showing weight increase over the 4 weeks


sum_data_Weight_W1$Week=1
sum_data_Weight_W2$Week=2
sum_data_Weight_W3$Week=3
sum_data_Weight_W4$Week=4
sum_data_weight_euth$Week=5

all_weight_sum <- rbind(sum_data_Weight_W1, sum_data_Weight_W2, sum_data_Weight_W3, sum_data_Weight_W4, sum_data_weight_euth)

library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:Hmisc':
## 
##     subplot
## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate, rename, summarise
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
w <- ggplot(all_weight_sum, aes(x = Week, y = mean, fill = drug)) +
  geom_bar(stat = "identity",position=position_dodge())+
  facet_wrap(~model, nrow = 1)+
                  labs(x = "", fill= "Intervention")

w+scale_fill_manual(values=c("#0000FF", "#FF0000"))

Difference in Weight (Euth - W1)

# weight difference 

full_data_nested$weight_diff <- full_data_nested$weight_euth - full_data_nested$Weight_W1
anova(lm(weight_diff ~ model * drug, full_data_nested))
## Analysis of Variance Table
## 
## Response: weight_diff
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## model       1 5343.2  5343.2 37.5706 1.293e-06 ***
## drug        1   61.9    61.9  0.4351    0.5149    
## model:drug  1  223.1   223.1  1.5690    0.2207    
## Residuals  28 3982.1   142.2                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_weight_diff <- ggplot(full_data_nested, aes(x=model, y=weight_diff, fill=drug)) + 
    geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8))
p_weight_diff + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), 
                  geom="errorbar", color="black", position=position_dodge(0.8), width=0.2 )+
    stat_summary(fun.y=mean, geom="point", color="black",position=position_dodge(0.8))+
  scale_fill_manual(values = c("#0000FF", "#FF0000"))+
    ylab("Weight (g) Gain across the study")+
  labs(x = "", fill= "Intervention")
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.